clc; clear all; close all
addpath(genpath('functions')) 
addpath('output_mat')


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~  1) Solve models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%% 1.1) Baseline model
clear all;
title = '01_baseline';  
baseline_parameters;   
steadystate; 
solve_01_baseline;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)
fprintf('Relative hours response (in %%):        %1.4f%% \n',mean((1-irf_mp_sss(:,pos.tfp)*100)./(1-irf_mp_sss(:,pos.n)*100)-1))

% Simulate
simulate_01_baseline;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos','irf_mp_simulated')


%% 1.2) Baseline model with alternative Taylor rule
clear all;
title = '01_alttaylor';  
baseline_parameters;     
sig_eps_i = 1.160*sig_eps_i;
steadystate; 
solve_01_alttaylor;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)


%% 1.3) Baseline model with high eta
clear all;
title = '01_highEta';  
baseline_parameters;     
eta       = 12; 
sig_eps_i = 0.00331; 
steadystate; 
solve_01_baseline;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)


%% 1.4) Baseline model with trend inflation
% 0.5% quarterly trend inflation
clear all;
title = '01_trend_infl_50'; 
baseline_parameters;     
pibar     = 1.0050; 
sig_eps_i = 0.00385; 
theta_1   = theta_2;
steadystate; 
solve_01_trend;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)

% 0.0% quarterly trend inflation
clear all;
title = '01_trend_infl_00'; 
baseline_parameters;     
pibar     = 1.0000; 
sig_eps_i = 0.00438; 
theta_1   = theta_2;
steadystate; 
solve_01_trend;
save(['output_mat/results_' title],'irf_mp_sss','irf_tfp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)


%% 1.5) Baseline model solved around deterministic steady state
clear all;
title = '01_hetDet';  
baseline_parameters;   
sig_eps_i = 0.0044;
steadystate; 
solve_01_detstst;
save(['output_mat/results_' title],'irf_mp_pos_sss','irf_mp_neg_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_pos_sss(pos.i_ann,:)))*100)


%% 1.6) Model with homogeneous price rigidity solved around deterministic steady state
clear all;
title = '01_homDet';  
baseline_parameters;  
theta_1   = theta;
theta_2   = theta;
theta_3   = theta;
theta_4   = theta;
theta_5   = theta;
sig_eps_i = 0.0028;
steadystate; 
solve_01_detstst;
save(['output_mat/results_' title],'irf_mp_pos_sss','irf_mp_neg_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_pos_sss(pos.i_ann,:)))*100)


%% 1.7) Model with specific labor
clear all;
title = '02_specLabor';  
baseline_parameters;  
sig_eps_i = 0.218*sig_eps_i;
steadystate; 
solve_02_specific;
save(['output_mat/results_' title],'irf_mp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)


%% 1.8) Model with Rotemberg price adjustment
clear all;
title = '03_Rotemberg';  
baseline_parameters;  
sig_eps_i = 0.00158;
fac     = 50;
theta_5 = 2;
theta_4 = fac*8;
theta_3 = fac*45;
theta_2 = fac*250;
theta_1 = fac*1950;
steadystate; 
solve_03_rotemberg;
save(['output_mat/results_' title],'irf_mp_sss','pos')

% Calibration targets
fprintf('Peak response of i (annualized, in %%): %1.4f%% \n',max(abs(irf_mp_sss(:,pos.i_ann)))*100)


%% 1.9) Model with menu costs or Calvo 
clear all;
title = '04_Menu';  
cost  = 1; 
baseline_parameters;  
discretization;
solve_04_menu_cost;
save(['output_mat/results_' title],'avg_probA','avg_dprice')

title = '04_Calvo';  
cost  = 2; 
baseline_parameters;  
discretization;
solve_04_menu_cost;
save(['output_mat/results_' title],'avg_probA','avg_dprice')






%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~  2) Create figures ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%% 2.1) Load all data

clc; clear all; close all
load results_01_baseline
irf_mp_baseline            = irf_mp_sss; 
irf_tfp_baseline           = irf_tfp_sss; 
irf_mp_baseline_sim        = irf_mp_simulated; 
pos_baseline               = pos; 
load results_01_alttaylor
irf_mp_baseline_alttaylor = irf_mp_sss; 
load results_01_highEta
irf_mp_highEta            = irf_mp_sss; 
load results_01_trend_infl_00
irf_mp_trend_infl_00      = irf_mp_sss; 
load results_01_trend_infl_50
irf_mp_trend_infl_50      = irf_mp_sss; 
irf_mp_trend_infl         = NaN([size(irf_mp_trend_infl_00),2]);
irf_mp_trend_infl(:,:,1)  = irf_mp_trend_infl_00;
irf_mp_trend_infl(:,:,2)  = irf_mp_trend_infl_50;
pos_trend_infl            = pos; 
load results_01_hetDet
irf_mp_hetDet_pos         = irf_mp_pos_sss'; 
irf_mp_hetDet_neg         = irf_mp_neg_sss'; 
pos_det                   = pos;
load results_01_homDet
irf_mp_homDet_pos         = irf_mp_pos_sss'; 
irf_mp_homDet_neg         = irf_mp_neg_sss'; 
load results_02_specLabor
irf_mp_specLabor          = irf_mp_sss; 
pos_specLabor             = pos; 
load results_03_Rotemberg
irf_mp_Rotemberg          = irf_mp_sss; 
pos_Rotemberg             = pos; 
load results_04_Menu
avg_probA_Menu            = avg_probA;
avg_dprice_Menu           = avg_dprice;
load results_04_Calvo
avg_probA_Calvo           = avg_probA;
avg_dprice_Calvo          = avg_dprice;
clear irf_mp_sss irf_tfp_sss pos avg_probA avg_dprice


%% 2.2) Create all figures

% Figure 5: a-f
plot_irf(irf_mp_baseline(:,pos_baseline.i_ann),[],'i',1,0,[-0. 0.3]);
plot_irf(irf_mp_baseline(:,pos_baseline.tfp),[],'tfp',1,0,[-0.4 0.0]);
plot_irf(irf_mp_baseline(:,pos_baseline.y),[],'y',1,0,[-.8 0]);
plot_irf(irf_mp_baseline(:,pos_baseline.agg_markup),[],'agg_markup',1,0,[-0.01 4]); 
plot_irf([irf_mp_baseline(:,pos_baseline.markup_1) irf_mp_baseline(:,pos_baseline.markup_2) irf_mp_baseline(:,pos_baseline.markup_3) irf_mp_baseline(:,pos_baseline.markup_4) irf_mp_baseline(:,pos_baseline.markup_5)],[],'markups',1,0,[-1 5]); 
plot_irf(irf_mp_baseline(:,pos_baseline.markupdisp),[],'markupdisp',0,0,[0 0.0015]);


% Figure 6: a-c
scale = irf_mp_baseline(1,pos_baseline.tfp)/irf_tfp_baseline(1,pos_baseline.tfp);
plot_irf((irf_mp_baseline_alttaylor(:,pos_baseline.y)-irf_mp_baseline(:,pos_baseline.y)),[],'alttaylor_gdp_diff',1,0,[-0.2 0]);
plot_irf_ax2(irf_mp_baseline(:,pos_baseline.markupdisp),scale*irf_tfp_baseline(:,pos_baseline.markupdisp),'tfp_markupdisp',0,1);
plot_irf(irf_mp_baseline(:,pos_baseline.agg_markup),irf_mp_highEta(:,pos_baseline.agg_markup),'agg_markup_eta',1,1,[-1 4]); 


% Figure F.1
plot_irf(irf_mp_baseline(:,pos_baseline.markupdisp),irf_mp_baseline_sim/100,'markupdisp_simulated',0,1,[0 0.0015]);


% Figure H.1
plot_irf_menu(avg_probA_Menu ,avg_dprice_Menu ,'precautionary_menu')
plot_irf_menu(avg_probA_Calvo,avg_dprice_Calvo,'precautionary_calvo')


% Figure H.2: a-f
plot_irf(irf_mp_baseline(:,pos_baseline.i_ann),irf_mp_baseline_alttaylor(:,pos_baseline.i_ann),'alttaylor_i',1,1,[-0. 0.3]);
plot_irf(irf_mp_baseline(:,pos_baseline.tfp),irf_mp_baseline_alttaylor(:,pos_baseline.tfp),'alttaylor_tfp',1,0,[-0.6 0.0]);
plot_irf(irf_mp_baseline(:,pos_baseline.y),irf_mp_baseline_alttaylor(:,pos_baseline.y),'alttaylor_y',1,0,[-1 0]);
plot_irf(irf_mp_baseline(:,pos_baseline.agg_markup),irf_mp_baseline_alttaylor(:,pos_baseline.agg_markup),'alttaylor_agg_markup',1,0,[0 4]);
plot_irf([irf_mp_baseline(:,pos_baseline.markup_1) irf_mp_baseline(:,pos_baseline.markup_2) irf_mp_baseline(:,pos_baseline.markup_3) irf_mp_baseline(:,pos_baseline.markup_4) irf_mp_baseline(:,pos_baseline.markup_5)],...
         [irf_mp_baseline_alttaylor(:,pos_baseline.markup_1) irf_mp_baseline_alttaylor(:,pos_baseline.markup_2) irf_mp_baseline_alttaylor(:,pos_baseline.markup_3) irf_mp_baseline_alttaylor(:,pos_baseline.markup_4) irf_mp_baseline_alttaylor(:,pos_baseline.markup_5)],'alttaylor_markups',1,0,[-1 6]);
plot_irf(irf_mp_baseline(:,pos_baseline.markupdisp),irf_mp_baseline_alttaylor(:,pos_baseline.markupdisp),'alttaylor_markupdisp',0,0,[-0.00 0.0015]);


% Figure H.3: a-f
scale = irf_mp_baseline(1,pos_baseline.tfp)/irf_tfp_baseline(1,pos_baseline.tfp);
plot_irf(irf_mp_baseline(:,pos_baseline.i_ann),scale*irf_tfp_baseline(:,pos_baseline.i_ann),'tfp_i',1,0,[-0. 0.3]);
plot_irf(irf_mp_baseline(:,pos_baseline.tfp),scale*irf_tfp_baseline(:,pos_baseline.tfp),'tfp_tfp',1,0,[-0.4 0.0]);
plot_irf(irf_mp_baseline(:,pos_baseline.y),scale*irf_tfp_baseline(:,pos_baseline.y),'tfp_y',1,0,[-.8 0]);
plot_irf(irf_mp_baseline(:,pos_baseline.agg_markup),scale*irf_tfp_baseline(:,pos_baseline.agg_markup),'tfp_agg_markup',1,0,[-0.5 3]);
plot_irf_ax2([irf_mp_baseline(:,pos_baseline.markup_1) irf_mp_baseline(:,pos_baseline.markup_2) irf_mp_baseline(:,pos_baseline.markup_3) irf_mp_baseline(:,pos_baseline.markup_4) irf_mp_baseline(:,pos_baseline.markup_5)],...
             scale*[irf_tfp_baseline(:,pos_baseline.markup_1) irf_tfp_baseline(:,pos_baseline.markup_2) irf_tfp_baseline(:,pos_baseline.markup_3) irf_tfp_baseline(:,pos_baseline.markup_4) irf_tfp_baseline(:,pos_baseline.markup_5)],'tfp_markups',1,0);
plot_irf_ax2(irf_mp_baseline(:,pos_baseline.markupdisp),scale*irf_tfp_baseline(:,pos_baseline.markupdisp),'tfp_markupdisp',0,1);


% Figure H.4: a-f
plot_irf(irf_mp_baseline(:,pos_baseline.i_ann),irf_mp_highEta(:,pos_baseline.i_ann),'i_eta',1,0,[-0. 0.3]);
plot_irf(irf_mp_baseline(:,pos_baseline.tfp),irf_mp_highEta(:,pos_baseline.tfp),'tfp_eta',1,0,[-.8 0.0]);
plot_irf(irf_mp_baseline(:,pos_baseline.y),irf_mp_highEta(:,pos_baseline.y),'y_eta',1,0,[-.8 0]);
plot_irf(irf_mp_baseline(:,pos_baseline.agg_markup),irf_mp_highEta(:,pos_baseline.agg_markup),'agg_markup_eta',1,1,[-1 4]); 
plot_irf([irf_mp_baseline(:,pos_baseline.markup_1) irf_mp_baseline(:,pos_baseline.markup_2) irf_mp_baseline(:,pos_baseline.markup_3) irf_mp_baseline(:,pos_baseline.markup_4) irf_mp_baseline(:,pos_baseline.markup_5)],...
         [irf_mp_highEta(:,pos_baseline.markup_1) irf_mp_highEta(:,pos_baseline.markup_2) irf_mp_highEta(:,pos_baseline.markup_3) irf_mp_highEta(:,pos_baseline.markup_4) irf_mp_highEta(:,pos_baseline.markup_5)],'markups_eta',1,0,[-1 5]); 
plot_irf(irf_mp_baseline(:,pos_baseline.mctilde),irf_mp_highEta(:,pos_baseline.mctilde),'mc_eta',1,0,[-4 1]); 


% Figure H.5: a-f
plot_irf(irf_mp_hetDet_pos(:,pos_baseline.i_ann),irf_mp_hetDet_neg(:,pos_baseline.i_ann),'hetDet_i',1,1,[-0.29 0.30]);
plot_irf(irf_mp_hetDet_pos(:,pos_baseline.y),irf_mp_hetDet_neg(:,pos_baseline.y),'hetDet_y',1,0);
plot_irf(irf_mp_hetDet_pos(:,pos_baseline.markupdisp),irf_mp_hetDet_neg(:,pos_baseline.markupdisp),'hetDet_markupdisp',0,0);
plot_irf(irf_mp_homDet_pos(:,pos_baseline.i_ann),irf_mp_homDet_neg(:,pos_baseline.i_ann),'homDet_i',1,1,[-0.29 0.30]);
plot_irf(irf_mp_homDet_pos(:,pos_baseline.y),irf_mp_homDet_neg(:,pos_baseline.y),'homDet_y',1,0);
plot_irf(irf_mp_homDet_pos(:,pos_baseline.markupdisp),irf_mp_homDet_neg(:,pos_baseline.markupdisp),'homDet_markupdisp',0,0);


% Figure I.1: a-f
plot_irf(irf_mp_specLabor(:,pos_specLabor.i_ann),[],'specLabor_i',1,0,[-0. 0.3]);
plot_irf(irf_mp_specLabor(:,pos_specLabor.tfp),[],'specLabor_tfp',1,0);
plot_irf(irf_mp_specLabor(:,pos_specLabor.y),[],'specLabor_y',1,0,[-.8 0]);
plot_irf(irf_mp_specLabor(:,pos_specLabor.agg_markup),[],'specLabor_agg_markup',1,0);
plot_irf([irf_mp_specLabor(:,pos_specLabor.markup_1) irf_mp_specLabor(:,pos_specLabor.markup_2) irf_mp_specLabor(:,pos_specLabor.markup_3) irf_mp_specLabor(:,pos_specLabor.markup_4) irf_mp_specLabor(:,pos_specLabor.markup_5)],[],'specLabor_markups',1,0);%,[-1 5]); 
plot_irf(irf_mp_specLabor(:,pos_specLabor.markupdisp),[],'specLabor_markupdisp',0,0);


% Figure I.2: a-f
plot_irf(irf_mp_Rotemberg(:,pos_Rotemberg.i_ann),[],'Rotemberg_i',1,0,[-0. 0.3]);
plot_irf(irf_mp_Rotemberg(:,pos_Rotemberg.tfp),[],'Rotemberg_tfp',1,0);
plot_irf(irf_mp_Rotemberg(:,pos_Rotemberg.y),[],'Rotemberg_y',1,0,[-.8 0]);
plot_irf(irf_mp_Rotemberg(:,pos_Rotemberg.agg_markup),[],'Rotemberg_agg_markup',1,0);
plot_irf([irf_mp_Rotemberg(:,pos_Rotemberg.markup_1) irf_mp_Rotemberg(:,pos_Rotemberg.markup_2) irf_mp_Rotemberg(:,pos_Rotemberg.markup_3) irf_mp_Rotemberg(:,pos_Rotemberg.markup_4) irf_mp_Rotemberg(:,pos_Rotemberg.markup_5)],[],'Rotemberg_markups',1,0);%,[-1 5]); 
plot_irf(irf_mp_Rotemberg(:,pos_Rotemberg.markupdisp),[],'Rotemberg_markupdisp',0,0,[0 0.0003]);


% Figure I.3: a-f
plot_irf(irf_mp_trend_infl_00(:,pos_trend_infl.i_ann),irf_mp_trend_infl_50(:,pos_trend_infl.i_ann),'trend_infl_i',1,1);
plot_irf(irf_mp_trend_infl_00(:,pos_trend_infl.tfp),irf_mp_trend_infl_50(:,pos_trend_infl.tfp),'trend_infl_tfp',1,0);
plot_irf(irf_mp_trend_infl_00(:,pos_trend_infl.y),irf_mp_trend_infl_50(:,pos_trend_infl.y),'trend_infl_y',1,0);
plot_irf(irf_mp_trend_infl_00(:,pos_trend_infl.agg_markup),irf_mp_trend_infl_50(:,pos_trend_infl.agg_markup),'trend_infl_agg_markup',1,0);
plot_irf([irf_mp_trend_infl_00(:,pos_trend_infl.markup_1) irf_mp_trend_infl_00(:,pos_trend_infl.markup_2) irf_mp_trend_infl_00(:,pos_trend_infl.markup_3) irf_mp_trend_infl_00(:,pos_trend_infl.markup_4) irf_mp_trend_infl_00(:,pos_trend_infl.markup_5)], ...
         [irf_mp_trend_infl_50(:,pos_trend_infl.markup_1) irf_mp_trend_infl_50(:,pos_trend_infl.markup_2) irf_mp_trend_infl_50(:,pos_trend_infl.markup_3) irf_mp_trend_infl_50(:,pos_trend_infl.markup_4) irf_mp_trend_infl_50(:,pos_trend_infl.markup_5)],'trend_infl_markups',1,0);
plot_irf(irf_mp_trend_infl_00(:,pos_trend_infl.markupdisp),irf_mp_trend_infl_50(:,pos_trend_infl.markupdisp),'trend_infl_markupdisp',0,0);%,[0 0.0003]);%,[0 0.0015]);

